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Evaluating  A7  D  A  for  Sparse  liatricas:  Analysis 

Amnon  Gonen 
Naval  Postgraduate  School 

ABSTRACT 

The  evaluation  of  the  matrix  product  A7  A  or  A7  D  A  .where  A  .  \ 

is  an  mxn  real  matrix  and  D  an  rrixm  diagonal  matrix,  is  a  funda-  v 

mental  operation  for  many  algorithms.  Vie  analyze  the  evaluation  j  /» t  ^ •  •  - 

of  A7  A  for  several  configurations  of  sparse  matrices  A  ail  of  which  j  ^  '  - 

have  the  same  sparsity.  The  complexity  of  the  evaluation  is  ;  ' 

i  . .  i,  • 

estimated,  and  application  to  certain  problems  of  optimization  are 
given.  ■ 

**  ! 

Key  ‘fords:  Sparce  Matrix,  Hessian  evaluation.  Optimization 


L  EITF.ODTJ  CTION 

Many  fundamental  algorithms  in  numerical  analysis  include  the 
evaluation  of  AT  A  or  AJ-D  A,  where  A  is  a  real  mxn  matrix  and  D  is  a 
diagonal  mm  matrix.  Examples  are  given  in  papers  cn  factorization  cf 
matrices  or  problems  of  minimization  in  which  the  Hessian  has  this  form 
(  Gay  [l]  Gotten  &  Avrie!  [3]  ).  The  extended  use  of  this  product  motivates 
the  question  of  reducing  its  complexity, 

This  research  ttm  pertially  supported  by  the  XPS  Foundation  Research  Program. 
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The  purposes  of  this  paper  are: 

1.  To  relate  the  computational  complexity  of  Ar  D  A  to  the  sparsity  rate 
of  the  matrix  A  . 

2.  For  a  given  sparsity  rate,  to  distinguish  between  the  worst  and  the 
best  case. 

3.  To  provide  an  application  of  these  results. 

The  problem  of  multiplying  a  transpose  of  a  sparse  matrix  by  itself 
was  discussed  in  several  books  and  papers  e.g.  George  3c  Liu  [2]  in  which 
they  include  the  number  of  operations  required  for  this  multiplication. 
Gustavson  [4]  proposed  an  optimal  algorithm  for  multiplying  two  sparse 
matrices  A  B  where  AzRnxm  and  3zRm*--,  proving  that  the  number  of 
multiplication  N  satisfies  0 5=  /V < nmk  .  However,  the  connection  between 
the  number  of  operation  and  the  sparsity  rate  of  the  matrices  was  not 
discussed. 

Apparently,  it  seems  that  this  question  has  only  theoretical  meaning 
since  the  matrix  A  is  provided  and  therefore  the  number  of  operations  is 
known.  However,  in  this  paper  we  will  see  there  exist  some  cases  in  which 
the  configuration  of  this  matrix  A  can  be  designed  by  the  user.  In  these 
cases  it  make  sense  to  analyze  this  product  in  order  to  reduce  the 
number  of  operations. 

In  section  2  of  this  paper,  we  present  the  computational  complexity 
of  Ar  D  A  for  several  sparsity  patterns  of  A.  In  this  section,  we  establish 
our  results  on  the  assumption  that  the  number  of  nonzero  elements  of 


the  matrix  A  is  provided.  We  demonstrate  the  best  and  the  worst 
showing  that  in  the  best  case,  the  nonzero  elements  are  divided  homo¬ 
geneously  among  the  rows  of  A,  while  in  the  worst  case,  these  nonzero 
elements  are  confined  in  a  limited  number  of  rows. 

In  section  3  we  provide  an  example  from  optimization  theory,  in 
which  the  matrix  A  is  dense  and  by  applying  the  results  of  section  2  —a 
minimize  the  number  of  multiplication  in  the  evaluation  of  the  Hessian. 

In  this  paper,  all  vector  spaces  are  finite  dimensional  and  vectors  are 
column  vectors.  The  space  of  ail  nxrz  matrices  is  denoted  by  ;  the 
nonnegative  crthant  of  the  Euclidean  space  !x~  is  denoted  by  ~v  ;  the 
subset  of  all  integer  vectors  in  Rn  is  denoted  by  and  its  nennezauve 
orthant  by  /?.  For  a  matrix  A  we  denote  by  rH.  and  the  i-th  row  and 
the  ;'-t’n  column  respectively.  The  transpose  of  A  is  denoted  by  AT .  By 
the  norm  Jjzjj  we  mean  the  Euclidean  norm.  For  a  real  number  r  its 
integer  part  is  denoted  by  [rj.  Finally,  the  number  of  elements  in  the  set 
B  is  denoted  by  ]B)  ,  and  the  number  of  zero  elements  in  a  matrix  .1  is 
denoted  by  Z(A). 

2.  THE  CQiJPUTATICilAL  COIIPLEXITY  C7  AT  D  A. 

Let  A  be  in  Rmxn  vhth  N  nonzero  elements.  The  ratio  —  is  called 

mn 

the  sparsity  rate  of  the  matrix  A  and  denoted  by  o(h).  In  this  section  we 
assume  that  the  sparsity  rate  of  the  matrix  AeR™**  is  provided  and  that 
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each  row  of  A  includes,  at  least,  one  nonzero  element.  We  concentrate  on 
the  sparsity  pattern  of  A,  looking  for  the  best  and  the  worst  cases,  by 
means  of  the  number  of  operations  required  to  compute  AT  D-A  where  D 
is  a  diagonal  matrix  D^RmXn .  We  begin  our  exploration  in  the  worst  case, 
in  which  the  configuration  of  A  implies  the  maximum  number  of  multipli¬ 
cation.  Let  us  denote  by  the  number  of  nonzero  elements  in  c^.  ,  thus 

2  mi=N.  (2.1) 

t=i 

Our  first  Lemma  provides  us  the  number  of  operations  (multiplica¬ 
tions)  required  to  accomplish  the  product  AT-A. 


Lemma  2.1:  Let  AzRmyn  be  a  given  sparse  matrix,  then  the  product  A7  A 
can  be  computed  using 

>SE”V("*t  +  l)  (2.2) 

i  =  t 


multiplications. 

Proof:  The  product  A7  A  can  be  rewritten  as  a  sum  of  m  rank  1  matrices 


ArA=  £  ^ 

i=  1 


CL1. 


The  rank  1  matrices  a^.  a7,  are  symmetric.  Each  nonzero  element  a^j  of 
the  vector  Oi»  is  multiplied  by  all  ether  nonzero  elements  a^.  for  ,fcs :j  . 
Therefore,  the  number  of  multiplication  is 


t«i 


(2.4) 


-  0  - 


combining  (2.3)  with  (2.4)  yields  Lhe  proof  of  the  lemma. 


From  the  proof  above,  it  can  easily  be  seen  that  the  number  of  addi¬ 
tions  are  approximately  the  same  as  the  number  of  multiplications  since 
each  term  a*.*  a.:.}  is  accumulated  into  the  result  matrix  C;  C  =  AT  A. 


Corollary  2. i:  Let  h€/?m*r-  be  a  sparse  matrix  and  a  diagonal 

matrix  then  the  product  A~  D  4  can  be  computed  by 


J +  1)  +  -V 

i  =  l 


multiplications. 

Prcci:  T.Ve  first  compute  A  =  D  A  which  requires  N  multiplications  and 
then  substituting  a£  by  aJ.  in  (2.3)  yields  the  proof  of  the  corollary. 

In  order  to  find  the  sparsity  partem  which  yields  the  worst  case,  we 
have  to  maximize  (2.2)  provided  (2.1)  and  all  rru  are  positive  integers. 
Since  the  difference  between  (2.2)  and  (2.5)  is  .V,  it  is  enough  to  explore 
the  worst  case  for  the  product  Ar  A  that  will  yield  the  same  result  for 
A7  D  A.  Consequently,  a  new  problem  can  be  formulated  as  follow. 


(Al) 


max. 


77V  (m*  +  1) 
2 


(2.6) 


subject  to  the  constraint 
2*n-t  =jV 

t«i 
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and 


sn  ;  mjS:/1  (2.7; 

m 

This  problem  can  be  reduced  to  maximizing  £  m?  under  the  same  con- 

t=i 

straints.  Defining  xt  =  m*  -1  yields  the  following  problem 

(A2)  maxi  I*  I  I2  (3.8) 

subject  to  the  constraint 

f]x i  =  N-m  (2.9) 

t=i 


and 


O^Xt  <n  - 1 .  x.  e/1 


(2.10) 


We  will  prove  that  since  the  objective  function  is  convex  its  maximum  is 
attained  at  a  boundary  point.  An  integer  vector  xc/n  is  called  a  boundary 

point  of  problem  (A2)  if  there  exists  a  set  /  =  j;, . =  . ml  and  a 

unique  L-J  such  that 


n-1  ie/ 

(N-7n)--$(n-  1)  i-j 0 
0  0  the-nriss 


(2.11) 


where 


N  -m 


n-1 


(2.12) 


In  this  case  the  vector  2=x+s  where  e=(l . 1)  is  a  boundary  point  of 

problem  (Al). 


Fortunately,  from  the  symmetric  property  of  the  objective  function, 
the  optimal  value  does  not  depend  on  the  selected  boundary  point. 


Hence 

§x,2  =  h(n-l)2  +  [(/V-m)-h(n-l)]2  (2.13) 

1=1 

To  prove  that  (2.11)  is  the  solution  of  problem  (.42)  we  need  the  fol¬ 
lowing  lemma 


Lemma  2.2  :  Consider  the  integer  problem 

(A3)  maxSjxji3  (2.11) 

xsla 

subject  to  the  constraints 

=  (2.15) 

i=l 

0<^<.V  (2.13' 

where  X  and  M  are  positive  integers,  !,1  <,  K.  Problem  (A3)  has  a  solution 
x *  satisfying 


\x*\  | z  =  i5.72  +  (K - h.t/)~ 


(2.  r 


K 

where  i?  is  the  integer  part  of  -^-denoted  by 


X 


,  if  and  only  if 


(2.13) 

Moreover,  if  (2.18)  holds  then  every  solution  of  problem  (A3)  satisfies 


Proof:  It  is  immediate  that  if  n-M  < K  then  there  is  no  feasible  solution 
to  problem  (A3).  Therefore,  let  us  assume  that  (2.18)  holds  and  prove 
this  lemma  t_  induction  on  the  dimension  of  x.  If  n  =  l,  then  from  (2.15) 
we  have  x  =  K  ^  M  .  If  K  <  M  then  =  0  and  if  K  =  M  then  =  1  .  In  both 
cases  (2.17)  is  satisfied.  Assuming  the  assertion  is  true  for  all  n, 

Let  us  denote 


F( m ,  M ,  K)  =  maze  \\  \  x  \  ‘ 2:  V  xj  =  K ;  0  <=  a*  £  M ;  x  zF'  j 


to  i  :'s 


hence 


0=  Tcax,r[F(m-l,M,K~z^)+x^:xnzril  (2.2C) 


Since  by  the  induction  assumption,  (2.  IT)  holds  for  m  -1 


F{mM,K)=  maxrj\V!r!2  +  (K-xm-i>y)2+xZ  :  = 


K  —  x_ 


(2.21; 


L-P  .0  —  ^  i  F.-J  V,rVlpT'<3  ~  I  *—  [Kfjn  y  pf 

K  by  \f  )  Consider  the  maximization  problem  (2.21)  in  two  cases: 


1.  0<xm  sp. 


In  this  case  h  =  i3  and  the  problem  is 


max  \ tf/J*  +  (K-#U)2-Z(K -$M)- xn  +  2x£ :  xn?/J  j  (2.22 


o' 


Substituting  K-VM  by  p  ,  yields  the  maximization  of 


tfltfz  +  p3-2px.n+2xz 


(2. 


on' 


subject  to  the  constraint  0 «zmsp  and  xme/}  .  The  maximum  of  (2.23)  is 


attained  at  xm  =p  or  xm  -  0,  and 


|  |zj  |2  =  i>4/2  +  p2  =  1 5.1/2  +  (K  —$H)Z. 


(2.24) 


2.  p<xm-£M  . 

In  this  case  is  =  v  - 1  and  the  problem  is 

max  \  ('J  - 1  )M2  +  (p  +  II)2 -Z [p  -r  M)zn  +  2z4 

psj  nn:j 


(2.25) 


using  the  same  arguments  as  in  the  first  case,  the  maximum  is  attained 
at  xm  =  ',1 ,  thus 


|z  i  j 2  =  (tf  - 1  )/J*  +  Co  +  M)2  ~  2 pM  =  tfd/2  +p2  =  u.:;2  +  (X -*:.!)* 


<  0  [T£n; 


In  both  cases  (2.17)  holds,  which  complete  cur  proof. 


Applying  Lemma  2.2  to  problems  (Al)  and  (A2)  yields  the  following 
conclusion. 

Corollary  2.3:  Every  x<zn  satisfying  (2.  It)  is  a  solution  to  problem  (A2). 
Proof:  Suppose  zefj  satisfies  (2.11),  which  mean  that  (2.13)  holds.  Sub¬ 
stituting  M  =n-l  and  K  =  N-m  in  Lemma  2.2  implies  that  (2.17)  and 
(2.13)  arc  the  same,  and  Lemma  2.2  implies  that  z  is  a  solution  of  prob¬ 
lem  (A2). 


A  solution  to  problem  (Al)  can  be  established  by  setting 
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TTLi  = 


n 

N  —  m  —  i)  (n  —  1)  +1 
1 


ieJ 

i=j0 

otherwise 


(2.27) 


where  satisfies  (2.12)  and  ^  is  a  set  of  indices  \  J\  and  j>'o  is  an  index 
not  in  /.  The  computational  complexity  of  the  product  A7-.',  for  the  worst 
case  is  established  by  substituting  (2.27)  into  (2.6) 

p-ws  -  lil'Sn2  +  (N  -77i  -i3(n  -1)  +  l)2  +  m  -  -J  - 1  -r//].  (2.28) 


It  can  be  seen  that  in  the  worst  case  some  of  the  va-ch  achieve  tne 
upper  bound  n  ,  the  others  are  zero  and  only  one  of  the  rr-.-th  is  some¬ 
where  between  0  and  n.  This  mean  that  the  matrix  .1  has  as  many  full 
rows  as  possible,  the  rest  of  the  rows  have  one  element,  and  one  row  con¬ 
tains  the  remaining  nonzero  elements  of  .V. 

In  the  next  Lemma  a  new  bound  fcr  the  computational  complexity  is 
presented  which  enable  us  to  relate  the  sparsity  rate  and  the  mathemati¬ 
cal  effort. 


Lemma  2.4  The  computational  complexity  of  the  product  A 7  A  can  be 
bounded  by 

^^n(N-m)^ZN).  '  (2.2S) 

Proof:  Let  us  denote  by 

<p(k)-{k  7ia4-[(jV  -m)  -fc(n-l)  +  l]8  + m  -k  - 1  +  N).  (2.30) 


The  first  assertion  is  that 


-  li  - 


v®( 


-  m 


71  -  1 


\  -  /  /V  -  m , 

— r~>- 

71  —  1 


(2.31; 


If  we  denote  by  <f  =  — — 
calculation  yields  that 


N  —m 


n  - 1 


,  then  0<<f<l  •  A  straightforward 


r(~ — =  (.V  —m)(n  +  l)  +  m  +  N. 


Hence 


y(  ^  -ri^-7^  ~  f)  =(^V-77i) -(71+1) +771 +.V-  (2.53) 

=  ( .V -7n  )■  (n  + 1 )  +  m.  +  ,V - [  _7^(re--l)-f-m  -KV -r?ia  -.-{ f(n  —  1)  -r  l)2-'- 1  ]  = 

=  f7i3-«-(n-l)  +  l)3-f+  l  =  f  (1-c)  (n-1)2 

since  osf<l  the  last  expression  is  nonnegative  which  prove  car  f.rrt 
assertion.  The  rest  of  the  proof  is  established  by  Llie  following: 


Muc  =  Vtr  ( 


X  -  m 


n  - 1 


)*M  ~ 7~)  =  K[n ( ,V  - m )  +  2 ..v ] . 

n  —  l 


(2.341 


As  we  can  see,  (2.29)  provides  us  an  elegant  bound  for  the  computa¬ 
tional  complexity  of  the  worst  case.  This  bound  is  a  good  approximation 

to  the  computational  complexity  when  — is  close  to  its  integer  part. 

The  difference  between  the  mathematical  effort  of  computing  ATA  in  the 
worst  case  and  this  bound  is  actually  provided  in  the  right  hand  side  of 
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(2.33)  and  it  is 


(2.35) 


where  £  is  the  fraction  part  of 


N  -  m 
n  - 1 


The  bound  in  (2.2S)  can  be  expressed  as  a  function  of  the  sparsity 
rate  by  using  the  definition 


r~r 

aU)=^~ 

mn 


( 


which  leads  to  the  following  equality 

-m)  +  2X]  =  yznm{G{A)(ri  +  2)-l).  (2.37) 

It  is  interesting  to  observe  the  connection  between  the  bound  in  (2.23) 
and  the  mathematical  effort  to  accomplish  ATA  without  using  sparsity 
method  which  is 


\$n-(n  +  l)-m.  (2.38) 

The  difference  between  (2.33)  and  (2.29)  can  be  established  by  expanding 
these  two  formulas  achieving 

(n  +  i)m  —  }t[(A’ -m)n  +  2.V]  =  }£(n  -*-2 ){m  n  -.V).  (2.39) 

Dividing  and  multiplying  the  right  hand  side  of  (2.39)  by  mn  yield  the  fol¬ 
lowing  expression  for  the  difference 

yfcn  +  2)  m  n  {\  -  a{A))  (2.40) 


where  a(A)  is  the  sparsity  rate  of  A. 


2.2  The  best  case 
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In  our  discussion,  we  call  the  case  in  which  we  need  the  minimum 
number  of  multiplication  to  produce  ATA  provided  that  there  are  N 
nonzero  elements  in  A  the  best  case  .  The  number  of  operations  in  the 
best  case  can  be  derived  by  minimizing 

^nv(m1  + 1) 

L  p  v---**/ 

i=i  c 

subject  to  (2.S)  and  (2.7).  Yfithout  the  integer  restriction  it  is  immediate 
that  since  the  objective  function  is  convex,  the  solution  will  be  the  arith¬ 
metic  mean,  that  is,  for  all  i  ,  m..  =  — .  The  restriction  that  ail  the  :~a  have 

m 

to  be  integral  yields  the  solution 


N_ 

m 


(2.42) 


where  Z,  =  i  1 . ,  JqL  and  \L~J  -N- 


m.  Consequently,  th; 


number  of  multiplication  in  the  best  case  is 


_  m  l)  _ ,  . 

hie  —  uj  o  /it 

t*l  a 


_v_ 

m 


■V 


+  l)(iLY-mj^-| 


(2.4'J'i 


In  order  to  present  the  magnitude  of  the  difference  between  the  worst 

and  the  best  case,  let  us  assume  that  —and  ~ ^  are  integers.  In  this 

m  n  -  1 


case  (2.29)  holds  with  equality  and 
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Mte  ~xh  +  > n). 

m 

Subtracting  i±bQ  from  ^  yield 


(2.44) 


Mue  -M*e  =J$(niV-n77i 


m 


(2.45) 


=  ^n(iV-m)(l-(7(4)). 


If  vre  take,  for  example,  A'=J£m{n  +  l)  the  difference  v/ill  be 


v.'hile  juts  =  m^n  *  1HZ--~r  .  That  means  that  for  large  n,  is  approxi¬ 


mately  50%  more  than  ,ui5 . 


3.  APPLICATION 


In  this  section  we  present  an  example  in  which  the  product  A7  D  A  is 
required  where  D  is  a  diagonal  matrix  and  the  pattern  of  A  can  be 
designed  in  order  to  reduce  the  computational  effort.  Since  we  are  dis¬ 
cussing  the  number  of  zeroes  in  matrices,  let  us  denote  by  Z\A)  the 
number  of  zero  elements  in  the  matrix  A.  Consider  the  problem  Intro¬ 
duced  by  Gay  [  l] 

(PI)  min  ?(x)  =  £p.;(rt(r))  (Cm) 

i  =  l 

where  t*:/?71  ->R  ,  pi-.R^R  and  m^n.  Very  often  r(~)  =  (r,(x) . r_(x))  is  a 

linear  function  of  s  ,  (see  for  example  Gonen  S:  Avriel  [3],  or  the  least 
square  problem  in  Gay  [l])  which  mean 

r(x)=A-x-b.  (3.2) 

In  this  case,  the  gradient  and  Hessian  of  $?  have  particularly  simple  forms 

Vp(x)  =  ATp'(r(z))  (3.3) 

?3?(x)  =  A7  D  A  (3.4) 

where 

P'(r(z))  =  0'i(r,(z)) . ,p'm(rir.(s))]  (2.5) 

and 

D  =  diag[p"l(rl(z)) . p’7n(7-ri(x))]  (3.6) 

is  the  diagonal  matrix  with  diagonal  elements  p,\(ri(x))  .  Since  we  have  a 
simple  analytic  presentation  of  the  gradient  and  Hessian  ,  it  is  reasonable 


to  consider  using  Newton  method  to  construct  a  sequence  of  iterates 
which,  under  reasonable  conditions,  converge  to  a  local  minimizer.  This 
mean  that  the  product  ATDA  will  be  used  each  iteration  and  very  often 
this  computation  is  the  most  expensive  part  of  the  algorithm.  The  main 
idea  is  to  accomplish  an  initial  preparation  step  by  factoring 

A  =  BQ  (3.7) 

where  QzRr-'An  is  nonsingular  and  3z3:n*r-  has  (n2-n)  zeroes  in  it 
( Z(3 )  -  n).  The  next  step  is  to  substitute  by  y  in  (3.7)  leading  to  the 
problem 

(P2)  min  <p(x)  =  2pi(n(x))  (5. 3) 

where 

r(y)  =  3y-o.  (3.C) 

To  establish  the  connection  between  the  two  problems,  let  us  introduce 
the  following  Lemma: 

Lemma  3.1:  A  point  x *  satisfies  sufficient  conditions  for  minimum  of 
problem  Pi  with  r(x)  defined  by  (3.2)  if  and  only  if  y*  -  Qz  satisfies 
sufficient  conditions  for  minimum  or  problem  P2. 

Proof:  The  sufficient  conditions  for  minimum  of  problem  Pi,  where  r(x) 
satisfies  (3.2),  are: 


At  WAx')  -  0 


(3.10) 
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zT  A7  V-tp(Aj:‘)Az  >  0  (  J.  2 1 , 

for  all  z* 0.  Since  A=BQ  where  Q  is  a  nonsingular  matrix  (3.10)  is 
equivalent  to 

BTV?(By‘)  =  0  (3.12) 

and  (3.11)  can  be  rewritten  as 

zTQ7BT-Vz<p(£y’)3Qz  >  0  (3.13) 

for  all  z*0.  Since  Qz  -  0  ii  and  only  if  z  =  0  our  proof  is  completed. 

D 


It  is  important  to  mention  that  from  Lemma  3.1  we  can  deduce  that 
if  A  is  a  nonsingular  square  matrix  then  it  is  enough  to  minimize  p(y)  and 
the  minimizer  x’  will  satisfy  x‘  = 

In  our  next  lemma  we  introduce  a  set  of  matrices  Ae/?:n*n  such  that 
for  every  factorization  of  a  matrix  in  this  set:  A  =  BQ  where  O  is  a  non¬ 
singular  matrix,  the  matrix  B  will  have  at  most  n 2  -  n  zeroes 
(Z(B)£nz-n).  Next  we  show  a  practical  method  of  factorizing  a  full 
ranked  matrix  which  achieve  at  least  n2  -  n  zeroes:  in  general  we  cannot 
expect  more  . 

Lemma  3.2:  Let  AeRm*n  where  m>n  be  a  full  rank  matrix.  Let  fT  =  [/!,-/] 
be  an  m  by  n+m  matrix.  If  any  set  of  m  columns  of  A  are  linearly 
independent  then  for  every  factorization  A  =  BQ  where  Qs..Rn*n  is  a  non¬ 
singular  matrix  and  Bs.Rm*n,  the  matrix  B  will  include  ,  at  least, 


i  o 

-  lu  - 


n(m  +  l)-n2  nonzero  elements,  (that  is,  Z(B)  £  n2-n  ). 

Proof:  Consider  the  factorization  AQ~l  =  B  which  can  be  written  as  n 
identical  linear  systems 

A-(trl)~  -I-Bi-  0  j- 1 . n  (3.14) 

The  coefficients  matrix  £=[A,-/]  has  rank  m  and  any  mxm  submalrix  cf 

7i  has  full  rank.  Let  us  denote  by  a-  the  vector  b}  in  .  First  we 

[b*J) 

claim  that  x  has  at  least  m  +  l  nonzero  elements.  Suppose  x  has  less  then 
37i  + 1  nonzero  elements  then  it  has  at  least  n  zero  elements.  Suppose 
z<  =2.- =  ■  ■  ■  =2,  =0  and  define  C£Rm*m  to  be  a  submatrix  cf  X  with 

‘1  *2  *7% 

columns  Xy  where  jViu  for  all  According  to  the  lemma’s  assump¬ 

tion,  C  is  nonsingular  and  therefore  the  only  solution  to  C  v  =  0  is  y  =0 
which  mean  Q ^  is  zero.  This  contradicts  our  assumption  that  Q  is  ncn- 
singular.  Therefore  the  matrices  Q  and  B  together  have  at  least  nm+n 
nonzero  elements.  If  we  assume  that  all  the  zeroes  arc  in  B,  wc  still 
remain  with  n(m-i-l)— vz  nonzero  elements  in  B. 

■ 

Comment:  Any  Vandermonde  matrix  satisfies  the  conditions  of 
lemma  3.2  therefore  there  arc  infinitely  many  examples:  of  matrices  for 

which  cne  cannot  expect  to  get  more  than  ns  -n  zeroes  in  B. 

Next  we  introduce  a  practical  method  to  factorize  a  full  ranked 


matrix  A  with,  at  least,  nz-n  zero  elements  in  B. 
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The  factorization 

Let  A€/?m!<n  be  a  full  rank  matrix  where  m>n  .  Then  we  can  write 


A  = 


Ax 

Az 


(3.15) 


Suppose  .4,  is  nonsingular  nxn  matrix.  In  this  case  -we  can  take 


B  = 


/ 

A?  A  £ 


Q  =  [A,] 


(3.13) 


and  there  are  r:2— n  zeroes  in  3 .  However,  this  iactorizalio: 
case  cf  section  1.  In  order  to  accomplish  a  better  factor! 
assume  that  m  >2 n  in  this  case  we  can  write  the  matrix  as  fc 


ration,  1st  as 


„  -  N 


n  • 


where  A1€J?nxn  is  a  nonsingular  matrix,  .?3ef?~:<n  and  Assume 

that  Aa-Af 1  can  be  factorized  into  L  V  where  L  and  ”  are  '.ewer  and  upper 
triangular  matrices  respectively. 


B  = 


IT1 

i 


[Az~  A  i  ■  •  L 

L 


Q  =  U  Ax 


(3.  IS) 


will  give  us  a  factorization  with  n2-n  zeroes  in  B  and  its  form  will  be 
closer  to  uniform  distribution  of  the  zero  elements  among  the  rows  of  the 
matrix. 

It  is  interesting  to  observe  cases  in  which  the  matrix  A  is  not  cf  full 
rank.  We  will  shovr  that  in  some  cases  it  is  possible  to  achieve  mere  zeroes 
than  the  full  rank  case  and  in  other  cases,  the  opposite  is  true. 
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Lemma  3.3:  Let  /lef?171*11  where  rank  (A)  <  n.  A  sufficient  condition  icr  A  to 
have  a  factorization  A  =  BQ  ,  where  QzRnXn  is  a  nonsingular  matrix  and 
Bz.Rm%n  such  that  Z(B)  >  n2  -  n  is  that 

rank  (A)  +  n  <  m  +■  1  (3. 19) 


Proof:  Suppose  that  rank  (A)  -  k  ,  1  <.  k  <  n  .  Without  loss  of  generality  we 
may  assume  that  the  first  k  columns  of  A  are  linearly  independent  and 


the  last  (n-k)  columns  are  linear  combinations  of  the  first 
us  write  A  -  [A1,J42]  where  AitR771'*  and  There 


k  columns.  Lei 
exists  a  matrix 


£’e/tSx'n  such  that  As 


A yE  .  The  matrix  Ai  can  be  factorized  :c 


Ai  =  BiQx  according  to  (3.16)  where  B^P."**  has  zero 
QiZR',:*:t  a  ncnsixigular  matrix.  Let  3zR™m  be  the  matrix 


vls^snls  ~ 
v'ltn  S i  in  its 


first  K  columns  end  zeroes  in  its  last  (n-fc)  columns  and  let 


Q \  E 

i°  * 


fq  pa' 
C u  : 


Since  Q i  is  nonsingular,  Q  is  nonsingular  and  ,4  =  BQ.  In  this  case  B  has 
at  least  kz  -  k  +  m(n  -  k)  zeroes.  Recall  that  the  number  of  zeroes  in  3  in 
the  full  rank  case  is  n2  -  n  ,  it  follows  that  fc2  -  fc  +  m(n  -  k)  >  n2  - n  iff 
kz  -  k(m  +  1)  +  n(m  +  1)  -  n2  >  0  iff  A:2  -  n2  >  (m  +  1  )(k  -  n)  Since  k  <  n 
the  last  inequality  will  hold  iff  A:  +  n  <  m  +  1.  This  inequality  is  the 
sufficient  condition  in  (3. 19). 


Conclusions 


We  have  seen  in  this  paper  a  class  of  optimization  problems  for  ..nun 
the  Hessian  matrix  can  be  written  as  A7  D- A  where  AeRm*n  and  DsR:r-'xrn  a 
diagonal  matrix.  Vie  showed  that  in  several  cases,  the  matrix  A  can  be 
partially  designed  by  the  user  in  order  to  reduce  the  number  of  nonzero 
elements  to  a  minimum.  In  previous  sections  we  explored  the  pattern  of  a 
sparse  matrix  with  a  given  number  of  nonzero  elements.  V.'e  shewed  that 
in  order  to  minimize  the  computational  complexity  of  A"-D  A  we  should 
divide  the  nonzero  elements  uniformly  among  the  rows  of  A  and  ::  the 
nonzero  elements  are  confined  in  certain  rows  then  the  computational 
complexity  is  maximized. 

The  difference  between  the  evaluation  of  the  product  A7-A  by  method 
of  dense  matrices  and  the  upper  bound  fer  the  worst  case  using  spars? 
method  is  presented  in  (2.40').  It  can  be  seen  that  this  difference  depend? 

linearly  on  the  proportion  of  zero  elements  in  the  matrix  which  is 

77171 

.  Furthermore,  the  saving  in  using  sparse  method  is,  at  least,  ;;(ti +2)rr.r. 
times  this  proportion.  Since  }£(n + 2) rr.n  and  (2  30)  are  both  close  for  large 
m  and  n,  the  saving  is  at  least  the  number  of  operations  for  the  dense 
case  times  the  proportion  of  the  zeroes  elements. 

Finally  we  demonstrated  a  practical  method  for  factorizing  a  full 
ranked  matrix  A£Rm*n  into  B  Q  where  B  has  at  least  n2  -n  zero  ele¬ 
ments.  Furthermore,  we  presented  a  class  of  matrices  A  for  which  you 
cannot  expect  to  get  more  than  n2  -n  zero  elements. 


-  ;.7j  _ 

Unfortunately  ,  this  factorisation  is  not  optimal  since  Lhe  nonzero 
elements  are  not  distributed  uniformly  among  the  rows  and  this  question 
is  still  without  an  answer.  Secondly,  we  proved  that  we  can  achieve  at 
least  nz  -  n  zero  elements  in  B  if  A  is  full  ranked  or  rank  (A)  +  n  <m  +  1  . 
We  did  not  prove  anything  for  matrices  which  are  not  full  rank  and  do  not 
satisfy  (3.19)  .  The  author  conjecture  is  that  the  theorem  may  apply  also 


for  this  case. 
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